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ABSTRACT 

Comparisons of the integrated thermal pressure support of gas against its gravitational 
potential energy lead to critical mass scales for gravitational instability such as the 
Jeans and the Bonnor-Ebert masses, which play an important role in analysis of many 
physical systems, including the heuristics of numerical simulations. In a strict theoret- 
ical sense, however, neither the Jeans nor the Bonnor-Ebert mass are meaningful when 
applied locally to substructure in a self-gravitating turbulent medium. For this reason, 
we investigate the local support by thermal pressure, turbulence, and magnetic fields 
against gravitational compression through an approach that is independent of these 
concepts. At the centre of our approach is the dynamical equation for the divergence 
of the velocity field. We carry out a statistical analysis of the source terms of the local 
compression rate (the negative time derivative of the divergence) for simulations of 
forced self-gravitating turbulence in periodic boxes with zero, weak, and moderately 
strong mean magnetic fields (measured by the averages of the magnetic and thermal 
pressures). We also consider the amplification of the magnetic field energy by shear 
and by compression. Thereby, we are able to demonstrate that the support against 
gravity is dominated by thermal pressure fiuct nations, although magnetic pressure also 
yields a significant contribution. The net effect of turbulence in the highly supersonic 
regime, however, is to enhance compression rather than supporting over dense gas even 
if the vorticity is very high. This is incommensurate with the support of the highly 
dynamical substructures in magneto-turbulent fiuids being determined by local virial 
equilibria of volume energies without surface stresses. 
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1 INTRODUCTION 

Self-gravitating turbulent fluids cover an enormous range of 
length scales, from molecular clouds all the way up to cos- 
mological scales. Traditionally, the linear perturbation anal- 
ysis of an extended homogeneous medium by [Jeans ( 1902 ) is 
applied to infer the stability against gravitational collapse. 
Apart from a geometrical factor, the resulting critical mass 
and length scales also apply to isolated systems ("clouds") 
in dynamical equilibrium. This is a consequence of the viral 
theorem. The critical mass of a spherical isothermal cloud 
that is solely supported by its thermal pressure is known as 
Bonnor-Ebert mass (|Ebert||1955i Bonnor 1956). As pointed 



out by Ball esteros-Paredes| ( 2006| ) , however, the generalized 
virial theorem also includes the integrated energies of non- 
thermal motions and magnetic fields, as well as contributions 
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from stresses at the boundaries of non-isolated systems (see 
also Lequeux|2005 ) . The former can be interpreted as turbu- 
lent and magnetic pressures, respectively. While these mean 
quantities characterize the global state, no conclusions about 
local effects can be drawn. Nevertheless, it is commonly as- 
sumed that Jeans-like scales can be applied locally if the 
density is scaled from the mean value to peak values in over- 
dense structures. However, in self- gravitating turbulent gas, 
such structures are highly non-linear and generally not in 
equilibrium. Consequently, an extrapolation to this regime is 
questionable on theoretical grounds. Chandrasekhar ( 1951| 
incorporated turbulence into the Jeans criterion by intro- 
ducing an effective pressure that depends on the internal 



velocity dispersion of density enhancements (see alsojBonaz^ 
zola et al.|1987 ). This heuristic approach was later put onto 



a firmer basis by means of renormalization group theory 
( Bonazzola et al^|1992 ). However, a fundamental limitation 
remains even with this method: The Jeans mass with an ef- 
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fective pressure in the fashion of Chandrasekhar can be de- 
rived only for length scales above the energy injection scale 
of turbulence. To a certain degree, this might be the case for 
cosmological structure formation, where turbulence is typi- 
cally driven by gravity on length scales smaller than the size 
of halos (e g. |Sur et al.||201Q| [Federrath et al.||2Qllb| [Turk 



et al.||2012 Latif et al.||2013 ). Overdense structures in star- 
forming clouds, on the other hand, collapse on scales below 
the integral scale of turbulence, which can be much larger 
than the size of the whole cloud (Elmegreen Scalo|[200^ 



Elmegreen Burkert|20TQl [Klessen Hennebelle|2010t ."Of 

course, it is still plausible that gravity becomes dominant 
on length scales comparable to the local, density-dependent 
Jeans length, but a strict derivation is not possible if gravi- 
tational instabilities develop on time scales smaller than the 
dynamical time scale of turbulence. 

In this article, we attempt to overcome these limitations 
by considering the fully non-linear dynamical equations of 
self-gravitating gas. This approach necessarily involves sta- 
tistical analysis. Apart from the adopted simulation scenario 
and methodology, however, the results we obtain are inde- 
pendent of any theoretical assumptions. The dynamics of gas 
compression is characterized by a partial differential equa- 
tion for the divergence of the velocity. By applying the Pois- 
son equation, the Laplacian of the gravitational potential 
can be substituted by a term that is proportional to the lo- 
cal gas density. This term is amenable to a comparison with 
the other sources of gas compression, which are derived from 
the local thermal, turbulent, and magnetic pressures of the 
gas. If these source terms are positive, they drive expansion 
of the gas or slow its collapse. Otherwise, they aid to its 
compression. In a collapsing region, the rate of change is 
mainly driven by the gravitational compression rate and, 
consequently, the divergence becomes ever more negative 
(Schmidt 2009). The relevance of turbulent support relative 
to thermal support of the gas in the intracluster medium 



was already investigated by Zhu et al. (2010) and lapichino 
et al. (20TT|. We elaborate on this approach by generaliz- 
ing the divergence equation to magnet ohydro dynamics and 
by analyzing numerical data from simulations of strongly 
self-gravitating turbulence in the interstellar medium. 

The production of dense clumps in the turbulent two- 
phase medium resulting from colliding flows were investi- 
gated, for example, by Banerjee e t al.| ( |2009| . Turbulence 
produced by homogeneous and isotropic forcing in periodic 
boxes, on the other hand, has been demonstrated to be 
particularly suitable for calculating statistics on molecular- 
cloud scales (e. g., |Kritsuk et al. 2006, 2007] |Schmidt et al 



2009 Kritsuk et al 



2009| [Federrath et al.|2010b[|2011at . To 



avoid a globally inhomogeneous flow structure, we thus uti- 
lize data from adaptive mesh refinement (AMR) simulations 
of self-gravitating hydrodynamical (HD) turbulence (Krit- 
|suk et al. 2011a) and magn etohydrodynamical (MHD) tur- 
bulence ( Collins et al.|2012| . These simulations were initial- 
ized by forced isotropic turbulence without self-gravity and 
then AMR was used to compute the contraction of the gas 
under the action of gravity. Since the statistically homoge- 
neous forcing completely randomizes the turbulent flow, also 
the dependence on initial conditions is minimized. Particu- 
larly the HD simulation offers an extremely large dynamical 
range, which resolves collapsing objects very well. An impor- 
tant question regarding the simulation of self-gravitating gas 



is the utilization of sink particles. If excess mass is dumped 
into sink particles, the gas density can be kept below a given 
threshold even at the maximum refinement level. Such a 
threshold is typically set to a multiple of the Jeans length 
( Truelove et al. 1 1 9*97 ) . By following this approach, Federrath 
et al. (2010a) demonstrated that gravitationally collapsing 



cores in a turbulent cloud can be identified with a set of 
heuristic rules, which test the divergence of the flow, the lo- 
cal gravitational potential, the local Jeans length, etc. How- 
ever, two uncertainties remain. Firstly, the removal of mass 
from grid cells inevitably produces small- scale perturba- 
tions in the solution of the equations of gas dynamics. This 
is potentially problematic for the calculations presented in 
this article because they depend on higher-order derivatives 
of gas-dynamical variables, especially, in the high-density 
regions that would produce sink particles. Secondly, exist- 
ing sink particle prescriptions do not absorb magnetic fields. 
While magnetic field lines are at least partially dragged into 
collapsing objects, the gas absorbed into sink particles is ef- 
fectively decoupled from the magnetic field. In this work, we 
avoid these uncertainties by letting the gas contract to arbi- 
trarily high densities (limited only by the dynamical range 
and the robustness of the numerical simulations), although 
this comes at the cost of eventually violating the resolution 
limit set by a Truelove- like criterion. As a consequence, re- 
sults for extremely high densities are tentative. 

This article is structured as follows. In Section |2] we 
define support against gravity by the source terms in the 
partial differential equation for the divergence of the ve- 
locity field. The simulation methods and parameters are 
briefly summarized in Section |3] In Sections [4] and [5] we 
carry out a statistical analysis of the different source terms 
for purely hydrodynamical as well as magnetohydrodynam- 
ical turbulence simulations with weak and strong magnetic 
fields. Thereby, we are able to infer turbulent, thermal, and 
magnetic support depending on various gas-dynamical prop- 
erties. In addition, turbulent dynamo action and the com- 
pressive magnetic field amplification is analyzed. In the last 
section, we summarize our results and discuss their implica- 
tions. 



2 THE LOCAL SUPPORT FUNCTION 

The momentum equation for a perfectly conducting ideal 
fluid subject to the gravitational potential (j) and the mag- 
netic field B can be written as 

^ipv) + V ■ {pv^v) = -VP+ -(J X B) -pV(/). (1) 

ot c 

In ideal magnetohydrodynamics, the current density J is 
related to the curl of the magnetic field via Ampere's law, 



V xB = ^J, 

c 



(2) 



and the magnetic field is given by the compressible induction 
equation 



DB 



{B • V)v - Bd, 



where the substantial time derivative is defined by 
D d 



Dt dt 



(3) 



(4) 
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and d = V • V is the divergence of the velocity. 

The gravitational field g — — V0 directly produces di- 
vergence but not vorticity = V x v (the rotation 
operator applied to g is identical to zero). Because of the 
non-linear turbulent interactions of turbulent velocity fluc- 
tuations, however, rotational motions (vortices) affect the 
divergence of the velocity. For magnetized fluids, the fol- 
lowing equation for the rate of change of the divergence d 
is obtained by applying the divergence operator to the mo- 

120101) and by 



mentum equation (Schmidt 



substituting Ampere's law 
Dc? 1.2 ie|2 

1 



t][2009 



Zhu et al. 



V • ( V • Tn 



(5) 



+ — Vp-(VP- V-Tn,) 



The vorticity of the flow, u? = V x which is related to the 
antisymmetric part of the velocity derivative, tends to in- 



crease the divergence. The rate of strain, \S\ 
where 



(^2Sij Si^ 



dvj 
dxj 



+ 



dvj 
dxj 



^l/2 



(6) 



is the symmetric part of the velocity derivative, has the op- 
posite effect. The Maxwell stresses of the magnetic fleld are 
given by 



4n 



BiBj 



(7) 



The last term on the right-hand side corresponds to the mag- 
netic pressure Pm = P^/Stt, while the former term specifles 
the anisotropic tension of magnetic fleld lines. 

In a gravitationally collapsing region, flux conservation 
squeezes the magnetic fleld lines as long as the ideal MHD 
approximation holds. The rate of change of the magnetic 
pressure is easily obtained from the induction equation ([3|: 



Dt 



where S^j = Sij — ^d6ij is the trace- free rate-of-strain ten- 
sor (the contribution from the antisymmetric part Wij — 
Vi,j — Sij of the velocity derivative Vij vanishes because 
BiBjWij = 0). The flrst term on the right-hand side re- 
sults from the action of the turbulent shear on the magnetic 
fleld. This term can be either positive or negative. Turbu- 
lent dynamo action amplifles the magnetic fleld. The rate 
of increase of the magnetic pressure due to gas contraction 
(d < 0), on the other hand, is —^B^d. 

For given boundary conditions, the gravitational poten- 
tial (j) is determined by the mass density of all gravitating 
matter (baryonic gas, stars, dark matter) through the Pois- 
son equation. For a compact mass distribution surrounded 
by vacuum, the Poisson equation takes the simple form 



(9) 



For periodic boundary conditions, on the other hand, the 
mean density (p) is conserved and the gas contraction must 
be zero for uniform mass density in the absence of any flow 
and magnetic flelds. Consistency with the boundary condi- 
tions thus requires 



- 47tG{p - po) 



(10) 



so that 



: for p = Po = (p) , 



0, and B = 0. 



By using the notation of iZhu et al. ( 2010 ), we can gener- 



ically write the Poisson equation in the form 

= 47iGpoS. 



(11) 



Here, po denotes a suitable reference density, and it is under- 
stood that the dimensionless density variation 6 is consistent 
with the source term of the Poisson equation for any given 
system. For example, S = p/po — 1 for periodic boundary 
conditions. 

We can now write an equation for the rate of compres- 
sion of advected fluid elements: 
_ Dd 
Dt 



4:7tGpoS - A. 



(12) 



Each of the terms in this equation has the dimension of 
inverse time squared. By applying the product rule and V • 
B = to the magnetohydrodynamical terms in Eq. ([5]), it 
follows that the local support of the gas against gravity is 
given by the function 



A: 



\S\' 



+ 



dxidxi 

dp 



p + 



B^ 



1 dBj dBj 
4:7T dxj dxi 



(13) 



p2 dxi 



dxi 



1 

47r 



dXn 



The contributions of turbulence, thermal pressure, and mag- 
netic flelds to the support of the gas are given by 

1 d^p 



At 



\S\' 



and 



A„ 



Athe 



1 

4:7rp 
+ 



p dxidx, 



+ 



dxidxi 



-B' 



1 dp dP 

p2 dxi dxi ' 

dBi dBj 

H 

dxj dxi 



(14) 
(15) 



1 dp 



47rp^ dxi 



d 

dxi 



-B' 



P, 



dB^ 
dxj 



(16) 



respectively. The net effect is to resist gravitational collapse 

if A = Atherm + Amagn + Amagn > 0. NcgatlvC SUppOrt (A < 

0) can result, for instance, from strong shock compression 
(large rate of strain). For isothermal gas, P — c^p and co = 
const. The thermal support is then given by 



-Co In p. 



(17) 



The classical Jeans length can be obtained from a linear 



perturbation analysis of rate of compression (Eq. 12) for 
small density perturbations ^ ^ 1. Let the gas initially be 
at rest and the magnetic fleld be zero. Linearization of the 
continuity and induction equations then implies 



-d' and 



0, 



dt dt 

where d' is the divergence of the flow and B' the magnetic 
fleld induced by the density perturbation. With the usual 
plane- wave ansatz 5 oc exp['icjt + • cc], where uj is the an- 
gular frequency and k the wavenumber, we have d' = —iuj3 
and a dispersion relation is obtained by linearizing Eq. ( 12 ) 
for the compression rate: 



~2 
UJ 



-47tGpo + k^cl. 



(18) 
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Here, Cs is the speed of sound for the unperturbed state. 
Hence, the perturbation is unstable if 



Table 1. Simulation parameters. 



k < kj 



1/2 



(19) 



which is the result found by Jeans (1902). The Jeans cri- 
terion corresponds to an equilibrium between the thermal 
pressure support for a single wave mode and the gravity 



term in Eq. ( 12 ) 



7 2 2 

kjc^ 



An ad hoc modification of the Jeans criterion for a uni- 
form turbulent velocity dispersion was introduced by |Chan-| 
[dra sekha? (11951 ) and later derived in a rigorous manner 
by Bonazzola et al. ( |1992 ). However, their theory is valid 
only if turbulence is produced on length scales smaller than 
the Jeans lengthy which excludes most applications in astro- 
physics. 



Equation ( 12 ), on the other hand applies to density vari- 
ations of any magnitude in magnetized turbulent gas. We 
can distinguish the following different regimes: 

(i) For A > AnGpod > 0, the gas is supported against 
gravity. Since Dd/Dt > 0, gas contraction {d < 0) is slowed 
down. 

(ii) Fluid elements undergo transient phases of vanishing 
support (A <^ 47vGpo6) or non- gravitational compression 
(A < 0), e. g., by shocks, as long as turbulent fluctuations of 
the pressure, the velocity, and the magnetic field push the 
gas back into regions with positive support. However, as the 
overdensity 6 grows, it becomes increasingly improbable that 
a fluid element can escape the pull of gravity and collapse 
may ensue. 

(iii) For a fluid element in a collapsing region, d < 0, 
47vGpo6 ^ |A|, and the free-fall time scale ^ {47vGpo6)~'^^'^ 
associated with the gravity term becomes the dynamically 
dominant time scale. 

(iv) At the end of collapse, the gas reaches an equilibrium 
state, in which the volume- averaged support is balanced by 
gravity, i. e., (A) ^ {47vGpo6). For an isolated object, this 
corresponds to the simple virial equilibrium (possibly, with 
non-thermal kinetic and magnetic contributions). 

Since the terms in the local support function A are second 
derivatives or products of derivatives, this function has a 
factor k^ in Fourier space, which is analogous to the above 
Jeans equilibrium condition. Generally, however, there is 
whole spectrum of perturbations, ranging from the small- 
est to largest wave numbers of the system. 

As an example, Fig.[l]shows slices of the functions (14), 
( 15 ), and ( 16 ) normalized by 4:7TGpo for the isothermal MHD 
simulation with /3o — 20 (see Collins et al.|2012 ). The slices 
are centered on the density maximum. Also shown are the 
thermal and the magnetic pressur es as well as the den stro- 
phy Qi/2 = ||V X {p^^^v)\^ (see [Kritsuk et al.||2007| . We 



chose fourth-order finite differences to compute gradients 
and LaplaciansQ Large values of the support function are 



^ By using a higher order than the spatial order of the numeri- 
cal scheme that is used to integrate the fluid dynamics, we keep 
the error in the computation of the gradients smaller than the 
intrinsic error of the numerical solution. 
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Notes: (1) simulation set; (2) solver used; (3) RMS Mach num- 
ber, A4; (4) maximum number of AMR levels, A^; (5) reflnement 
factor, R] (6) minimum number of zones per Jeans length, Lj; 
(7) ratio of compressive to solenoidal motions, x'l (8) was driving 
continued through the collapse phase. 

clearly associated with steep pressure and denstrophy gra- 
dients. For isothermal turbulence, the thermal pressure is 
proportional to the density, while the denstrophy indicates 
both density fluctuations and vortices. While negative tur- 
bulent support is most prominent, the thermal and magnetic 
support swap their signs across shock front. Even in the cen- 
tral, refined region, where the gas density reaches its peak, 
the support is not predominately positive. These qualita- 
tive observations are confirmed by the statistical analysis in 
Section O 



3 SIMULATIONS 

The hydrodynamic simulations we analyse here were earlier 
presented in Padoan et al. (2005) and Kritsuk et al. ( 2011a|), 
and the MHD simulations were presented in [Collins et aL] 
(2012). Both sets of simulations were run with the Enzo 
cod^EI Table [l] summarizes the differences between the two 
simulations. Both simulations began with 512^ grids, with 



uniform density. Driving was done as in Mac Low (1999), 
wherein small static velocity perturbations are added every 
time step such that kinetic energy input rate is constant. 
Once a steady state was reached, gravity and AMR were 
initiated. In both simulations the power in the driving field 
is distributed approximately isotropically and uniformly in 
the interval k G [1,2]. 

The hydrodynamic simulation used the third order 
piecewise parabolic method ( Colella Woodward||1984 ). It 
was driven with a mixture of 60% solenoidal and 40% com- 
pressive motions, and the energy injection was such that 
the RMS Mach number was approximately 6. The inertial 
range ratio of compressive to solenoidal motions, and its de- 
pendence on Mach number and magnetic fields, is discussed 
in [Kritsuk et al.| ( [20To| . At t = 0, when gravity and AMR 
the driving is halted, though the kinetic 



are switched on, 
energy decay in the short duration of the collapse phase is 
insignificant. Five levels of AMR were used, with a refine- 
ment ratio of 4, and a refinement criterion such that the 
Jeans length was resolved by at least 4 zones. 

The MHD simulations used the second order piecewise 
linear method (Li et al. 2008), the divergence- free con- 



strained transport method described in [Gardiner &: Stone] 
(I2005D, and the divergence- free interpolation method of 



Balsara (2001). The code is described in detail in Collins 



Balsara 


( 


2001 


et al. (2010). ' 



The driving force was purely solenoidal, and 
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(e) Amagn/47rGpo (f) Pm = /Stt 



Figure 1. Slices of the dimensionless turbulent, thermal, and magnetic support (left column) and associated quantities (right column) 
for an MHD turbulence simulation with /5o = 20 at t = O.Stff. The size of the shown region is 0.3 times the box size. Black lines indicate 
the boundaries of refinement levels. 
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an RMS Mach number of approximately 9 was reached. 
Driving continued through the collapse phase of the simula- 
tion. Four levels of refinement were used, with a refinement 
ratio of 2, and refinement was such that the Jeans length 
was resolved by at least 16 zones. 



4 STATISTICAL ANALYSIS OF 

HYDRODYNAMIC TURBULENCE 

We use data from the hydrodynamical turbulence simula- 
tion performed by Kritsuk et al. (2011a) to calculate statis- 
tics of the thermal and turbulent support functions. The 
virial parameter of this model is a — 0.25 and the sonic 



Mach number is about 6. Unless otherwise stated, we use 
the final state of the simulation at time t — 0.43tff, where 
tff = (37r/32G'po)^^^ is the free-fall time for the mean density 

On the basis of the Jeans criterion, we expect that high- 
pressure regions are correlated with large, positive thermal 
support Atherm (Eq. 15). A quantity with the physical di- 
mension of pressure can be obtained by multiplication of 
Atherm with the gas density p and normalization with the 
square of the local grid scale A. Fig. |A1| in the appendix 



shows phase plots of the resulting quantity A pAtherm vs. 
P, i. e., two-dimensional histograms of the occupied volume 
fractions. Since the support functions are strongly fluctuat- 
ing quantities, logarithmic scaling is necessary to investigate 
the statistics over the full dynamical range of pressure or 
density. Consequently, we separate the support into positive 
and negative components A^pAthermi, where generically 



A4 



A if A ^ 0, 
otherwise. 



A_ 



-A 







if A ^ 0, 
otherwise. 



The phase plots of the thermal support show that both pos- 
itive and negative contributions are significant, but large 
values of Atherm + appear to be more frequent in comparison 
to Atherm-, particular at high pressures. 

For a quantitative comparison that allows us to discern 
trends, we compute profiles, i. e., mean values for small bins 
of pressure or other quantities (throughout this article, we 
use 0.05 dex as bin width; see Appendix [A] for further de- 
tails) . The mean positive and negative components add up to 
the mean net support, i. e., (/A) = (/A+) — (/A_) both for 
volume- and mass- weighted averaging and an arbitrary posi- 
tive function /. For A^pA, we calculate volume-weighted av- 
erages. Owing to the factor / = A^p, this effectively results 
in weighing the support by mass. This allows us to specify 
the typical magnitude of the support at a given pressure. In 
the top panel of Fig. |2] one can see that indeed the positive 
component of the thermal support dominates for high pres- 
sure. Since gravitational contraction produces much higher 
densities than the initial supersonic turbulence, the range 
of pressures also increases greatly. This can be clearly seen 
by comparing to the early instant t — O.ltff, for which the 
thermal support is also plotted. For t — 0.43tff, we can ap- 
proximate the pressure- averaged values by the asymptotic 
relation 

(A'Mtherm)p for P > 100. 

A value of unity corresponds to the thermal pressure at the 
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Figure 2. Averages of the positive and negative thermal sup- 
port A^pAtherm ± cis functions of the thermal pressure (top) and 
the averaged turbulent support A^pAturb± vs. enstrophy den- 
sity (middle) and denstrophy (bottom) at time t = 0.43tff. The 
positive and negative components are respectively shown as solid 
and dashed lines and the identity functions are indicated by a 
dot-dashed line in each plot. Factors of A^ are applied to obtain 
quantities with the dimension of pressure. For comparison, the 
light-coloured lines show the support functions for t = 0.1%. 



mean density. In this sense, a higher pressure corresponds 
to an enhanced support against gravity, as expressed by the 
formula for the Jeans length. But even for very high pres- 
sures, there is a non- vanishing fraction of negative support. 
Since P CoP, the above relation implies the typical mag- 
nitude Atherm (cq/A)^ at high densities. Naturally, the 
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Figure 3. Power spectrum of the turbulent support, computed 
with finite differences and pseudo-spectral derivatives (solid lines; 
the finite-difference computation is shown in red in the online 
version) . Also shown are the power spectra of the vorticity ou and 
the rate of strain. 



support functions depend on the grid resolution A. Never- 
theless, we can draw meaningful conclusions, as will become 
clear in the course of our analysis. For pressures or densities 
higher than about 10^, there appears to be a systematic de- 
viation of the thermal support from the above asymptotic 
relation. This becomes even more prominent in the statistics 
of the thermal support relative to gravity, as is shown below. 
The strong fluctuations for very high values of the pressure 
or density are due to the small samples in each bin. 

Let us now consider the turbulent support. A common 
interpretation of turbulent support is that turbulent veloc- 
ity fluctuations or, more specifically, eddy-like motions act 
against gravity. Let us assume the scaling laM[^^t'(^) ^ i'^^'^ 
for supersonic turbulent velocity fluctuations on the length 
scale i. Then we have (jj{£) ^ 6v{i)/i ^ for the mag- 

nitude of the corresponding vorticity. This implies that the 
dominant contribution to the numerically resolved vorticity 
comes from eddies at the (spatially varying) cutoff scale A, 
i. e., (jj oc A~^/^. Contributions from larger eddies are sup- 
pressed by ^ (A/^)^/^. The same argumentation applies to 



the rate of strain \S\ and it follows from Eq. (14) that the 
local turbulent support Aturb A~^. Since the turbulent 
pressure has the dimension of density times squared veloc- 
ity fluctuations, we consider pA^ Aturb, where the first term, 
lA^pcj^ = A^Q, can be interpreted as the turbulent pres- 
sure due to eddies on the grid scale and the second term, 
— IA^pIS*!^, corresponds to a negative pressure that is pro- 
duced by the strain. The quantity Q — \p(jO^ is the enstrophy 
density. Since the trace of the rate-of-strain tensor is the di- 
vergence large strain is particularly caused by shocks. For 
the gravitational support, the crucial question is whether 
the positive or the negative component is dominant. 



Profiles of the turbulent support. 



A pAturb- 



A>At. 








if > \S\, 
otherwise, 

if < \Sl 
otherwise. 



are plotted in Fig. [2] Regardless of the average enstrophy 
density, we find that the negative component A^pAturb - is 
dominant. For a wide range of enstrophies, both components 
of the turbulent support follow very closely linear relations 
and the total turbulent support function is approximately 
given by the negative turbulent pressure: 

(AVAturb)A2Q ^ -A'Q if A'Q > L (20) 

The factor A^ in this relation does not trivially cancel out 
because values from different refinement levels with vary- 
ing A contribute to the average for a particular value of 
A^Q. The above relation also holds for the earlier instant 
t — O.ltff, except that the maximum vorticity is much 
lower. The above relation also implies that the typical mag- 
nitude of p\S\^ is linearly related to Q. For dimensional 
reasons, the turbulent pressure may alternatively be esti- 
mated by A^Q i/2; where O1/2 is the denstrophy (see Sect. [2]). 
As argued by Kritsuk et al. (2007), the denstrophy com- 
bines the effects of eddies (through the rotation of v) and 
shocks (through the density gradient). The resulting pro- 
file in Fig. [2] (bottom panel) shows a very similar trend and 
also a nearly linear asymptote. The phase plots in Fig. |A2| 
show a substantial overlap between positive and negative 
turbulent support so that locally the support by vortices 
can be stronger than the compressive effect of turbulence. 
For a given vorticity, however, (A^pAturb +)a2q could ex- 
ceed (A^pAturb -)a2q only if the rate of strain were smaller 
than the vorticity in the majority of cells. 

The length scales, which are dominated either by strain 
or vorticity, can be inferred by decomposing the turbulent 
support in Fourier space. The scaling law uj{t) ^ im- 
plies a flat power spectrum of Aturb Q Since the turbulent 
support is quadratic in the velocity derivative, 

1 , 



Af 



\S\') 



turb 



(fc) 



kikj{viVj -\- ViVj)k dQfc 



the computation a power spectrum is straight-forward: 
1 

"2 

1 



|fc|=fe 



+ Vi. 



Sij Sij 



(21) 



In the above expressions, the Fourier transform of a field / 
is denoted by / and /* is its complex conjugate. The inte- 
grands are symmetrized by complex conjugation and inte- 
grated over spherical shells of radius k in Fourier space. By 



^ In code units with appropriate normalization. The actual ex- 
ponent of the scaling law is not important for our reasoning. It 
can be anything between the Kolmogorov and the Burgers scaling 
exponents. 



^ The scale-dependence Aturb ~ A"-*^ corresponds to a linear re- 
lation with the wavenumber k, which cancels with the factor 
for a spectrum that measures the power per unit wave number. 
The turbulent velocity scaling 6v(i) ~ i-^^"^, on the other hand, 
corresponds to the spectrum E(k) oc k~'^. 
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Figure 4. Averages of the thermal (top) and turbulent (bottom) support relative to the gravitational compression rate in Eq. \12\ as 
functions of the overdensity {t = 0.43tff). As in Fig. [2] positive and negative components are shown as solid and dashed lines and the 
support functions at time t = O.ltff are plotted as light-coloured lines. The left and right columns of plots show the support functions 
for the full AMR data and the root-grid representations, respectively. 



using the root-grid representations of the data[^ we com- 
puted Pturb(^) both from the Fourier transforms of the 
finite- difference approximations to Vij and from the pseudo- 
spectral derivatives ikjVi. The results are shown in Fig. |3] 
The agreement between the two methods of computation 
is very good. Only toward the cutoff wavenumber, there 
is a deviation due to Gibbs phenomenon. Remarkably, it 
turns out that Pturh{k) < for all wave numbers, i.e., the 
compression produced by shocks overcompensates the sup- 
port by turbulent vortices. This can also be seen from the 
evaluation of the power spectra of lu and \S\. We find that 
Puj{k) is smaller than ^P\s\{k) for all wave numbers. More- 
over, Pturb(^) is nearly flat in between the forcing range 
at the smallest wavenumbers and the numerical dissipation 
range, as expected from the aforementioned scaling argu- 
ments. Apart from the damping close to the cutoff wavenum- 
ber, one can also see the typical bottleneck bump (see Krit- 
suk et al.|[2007p 



To investigate the magnitude of the local support func- 
tion A relative to the gravitational compression rate 47vGpo6 



1000 



Atherm+/4;rGpo^ 

Aturb-/47rGpo(5 



100 




100 200 



Figure 5. Plots of the positive (thick solid) and negative (thick 
dashed) components of the total support at O.ltff . Also shown are 
the dominant components of the thermal and turbulent support 
functions (thin lines). 



^ A consistent compuation of the spectra for the higher refine- 
ment level is not possible because of the missing high-wave num- 
ber modes of turbulence in coarser regions, which are not refined 
by the Truelove-like criterion. 



in Eq. (12), we compute averages of the ratio A±/47vGpo5 
for logarithmic bins of the overdensity 1 + ^ = p/po- For this 
ratio, we use the mass per cell as weighing function. First, we 
consider the thermal and turbulent support functions sepa- 
rately (see Fig. |4|. In the early phase {t = O.ltff), the the 
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1+^ 




1+5 

Figure 6. Plots of the same functions as in Fig. [s] both for the 
full AMR data (top) and the root-grid data (bottom) at time 
t = 0.43%. 



support functions simply decrease with density. The total 
support is plotted together with the dominant thermal and 
turbulent components in Fig. [5] The critical value of unity is 
approached for the highest densities (6 ^ 100), which indi- 
cates that the gas is self- gravitating. Since Aturh - /^TvGpod 
is large, however, the gas is mainly compressed by shocks. 
For t = 0.43tff, the statistics plotted in Figs. [4] and |6] clearly 
show that the positive component of the thermal support 
and the negative component of the turbulent support con- 
tinue to dominate for all ^ > 0. Compared to t = O.ltff, 
the support remains nearly unaltered for the lower range 
of densities, but an intermediate range 100 < (5 < 10^ has 
formed, in which the ratio of support to gravity is nearly 
constant and Atherm + is roughly balanced by Aturb - (see 
top panel in Fig.|6|. For this range of densities, the simula- 
tions resolve the physics of self-gravitating turbulence very 
well. In the phase plots of Atherm /47rG'po^ and Aturb /47rG'po^ 
(see Figs. A3 and A4), one can see that the strong ther- 
mal support at high densities is associated with rare strong 
fluctuations, while the negative turbulent support extends 
toward higher values than the positive component. At all 
densities, there is a significant volume fraction for which 
the support relative to the compression rate is small. This 
is simply a consequence of the support function fluctuating 
between positive and negative values. The gas is collapsing 
only if the mean support is sufficiently small. 

In Fig. [7| the positive and negative components of the 
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Figure 7. Linear plot of the average total support at time t = 
0.43tff , showing the positive and negative components as well as 
the net support. The grey shaded region indicates the region for 
which A < 4:7tGpo6, i. e., negative compression rate. The full 
AMR data are plotted on the left, the root-grid representation on 
the right. 



total support are plotted in linear scaling together with the 
profile of the net support A/ AnGp^d^ where A = A+ — A_. 
Except for the strong spikes, which probably correspond to 
collapsing filaments that start to feel strong pressure sup- 
port, the net support is negative and of the order of the 
gravity term. For the floor, we have |A| ^ AnGpod and a 
median value A/AnGpoS ^ —0.67. This indicates an accel- 
erated contraction of the gas, as supersonic gas compression 
triggers gravitational collapse. The absence of a dominant 
positive thermal pressure component in this range of densi- 
ties also supports the explanation for the deviation of the 
power-law index of the density PDF tail from —1.5 given in 
Kritsuk et al.| (|2011a). There, the pressureless collapse solu- 
(1969) is used to explain a steeper, ^ —1.7, 

At densities 



Penston 



tion of 

slope of the tail observed in the simulation, 
around 10^, however, the thermal pressure of the gas tends 
to overcompensate gravity. This does not necessarily mean 
that the gas is expanding (d > 0), but that the contraction 



slows down (D d/Dt > 0) toward higher densities. Kritsuk 



et al. (2011a) show that an overdensity of 10^ marks the 
transition to rotationally supported cores, which entails a 
decreasing importance of thermal support. This transition 
can be seen as a flattening of the power-law slope of the 
density PDF (see Fig. 1 in their paper). For S > 10^, both 
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the thermal and the turbulent support decrease relative to 
gravity. However, the gas mass per cell ceases to be a small 
fraction of the local Jeans mass for extremely dense gas (the 
Jeans length is comparable to the cells size at the highest 
refinement level for 5 ^ xlO^^ and, consequently, the gas 
dynamics is not well resolved |j 

A different picture emerges if we compute the support 
functions from the root-grid representation of the data. This 
amounts to a smoothing or filtering operation that averages 
the data at the higher refinement levels down to the much 
coarser cells at the root grid. The averaging over finer cells 
conserves mass and momentum. Compared to the highest re- 
finement level, the ratio of the linear grid scales is 4^ = 1024. 
Consequently, the smoothing of the densest structures in 
the regions with maximum refinement is substantial. This 
is refiected in the much smaller density range of the sup- 
port functions plotted on the right-hand side of Fig . [4) 
where we computed the derivatives in Eqs. (15) and 
from the smoothed density, velocity, etc. fields. Most im- 
portantly, we find that the K±/ AnGp^S gradually decreases 
to values below unity for both the thermal and the turbu- 
lent support, thus indicating that the densest gas is strongly 
self-gravitating. The transition occurs at densities around 
5 ^ 10^. Moreover, the net support is positive, as one can 
see in the plot with linear scaling shown in Fig. [t] (bot- 
tom plot). Consequently, the properties of the smoothed 
structures are markedly different from the collapsing struc- 
tures with negative net support that appear at higher refine- 
ment levels. Apart from that, the positive turbulent support 
shown in Fig. [4] increases relative to the negative component 
for 5 > 1000, which suggests that the smoothed-out effect 
of vorticity in the dense, self- gravitating gas becomes more 
important. This does not contradict the power spectra in 
Fig. [3] where the vorticity and the rate of strain are aver- 
aged for given wavenumbers. The volume weighted power 
spectra are not sensitive to the high density, low volume 
fraction self- gravitating structures. Here, we find a relatively 
large vorticity compared to the rate of strain for very high 
densities. A composite plot of the support functions for the 
root-grid data is shown in the bottom panel of Fig. [6] The 
large values of K±/ Ai^Gpod for 5 > 10^, which can be seen 
for the full AMR data in the top panel, are absent. The very 
large peaks of the thermal support at the higher refinement 
levels are likely due to strongly localized and to some degree 
under-resolved structures that tend to be over-pressurized. 
In contrast, gravity in refined regions is more dominant for 
a given density at the root grid because the mass elements 
are larger, whereas local fiuctuations tend to be averaged 
out. Thus, we observe the simple trend of gravity becoming 
increasingly dominant for increasing density. 



5 STATISTICAL ANALYSIS OF 

MAGNETOHYDRODYNAMICAL 
TURBULENCE 

To analyze the support in magneto-turbulent fiuids, we use 
data from the simulations by Collins et al. (2012| for a weak 



Table 2. Ratio of mean thermal to magnetic pressures for the 
MHD simulations. 





/3o 
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0.1 


20 


0.2037 


0.5 


20 


0.2056 


0.1 


0.2 


0.0336 


0.5 


0.2 


0.0339 



^ Regardless of the significance of the local Jeans length for 
gravitational instability, this is suggested by simple time-scale 
arguments. 



magnetic field with a mean thermal-to-magnetic pressure 
ratio = 20 and a moderately strong field with /3o = 0.2. 
The case of hydro dynamical turbulence considered in the 
previous section corresponds to the limit ^ 00. Since tha 
parameter fio specifies the strength of the initial, uniform 
field, we calculated the ratio of the mean thermal and mag- 
netic pressures, which are summarized in Table [2] As one 
can see, the magnetic pressure of the turbulent gas is much 
higher than the initial value. For brevity, we use the param- 
eter to tag the two MHD models. The virial parameter of 
both models is 1 and the sonic Mach number M is about 9. 
Since — v^^^/cl and /3 — 87Tpo/cs{B'^) , it follows that 
mean magnetic energy is comparable to the mean kinetic 
energy for /3o = 0.2 (corresponding to /3 ~ 0.034), while it is 
significantly smaller for /3o = 20 (corresponding to /3 ^ 0.2). 
Compared to the hydrodynamical model, gravity is weaker 
in the MHD simulations, but the turbulence energy is higher. 
Moreover, in contrast to the hydrodynamical simulation, we 
have a significantly lower limit of ^ ~ 6.35 x 10^ for resolved 
overdensities (see Sect.|3|. 

Since the impact of magnetic fields depends on the am- 
plification by shear and compressions, we first consider the 
two corresponding source terms in Eq. ^ for the time evolu- 
tion of the magnetic pressure. The first term arises from pure 
(trace- free) shear, while the second term is due to gravity 
and shocks (both are associated with highly negative diver- 
gence, which results in positive —B^d). The amplification of 
the magnetic field induced by the shear of the turbulent fiow 
can be understood as small-scale dynamo action. The mean 
positive and negative components of these terms are plotted 
as functions of the overdensity in Fig. [S] There is clearly a 
net amplification of the magnetic field in the high-density 
gas for both /3o = 0.2 and 20. However, in the initial phase 
(t = O.ltff), the total amplification is small because the am- 
plification at high densities (solid light lines) is roughly coun- 
terbalanced by the weakening of the field at low densities 
(dashed light lines). This indicates that the magnetic field 
is close to saturation, which makes sense because we start 
with statistically stationary turbulence at time t = 0. But 
the graphs for t — 0.5tff show that gravitational contraction 
causes further amplification in overdense gas (^ > 1), which 
is nearly matched by shear-induced amplification. In com- 
parison to /3o = 0.2, the field amplification in the weak-^e\d 
case (/3o = 20) is noticeably larger for overdensities in the 
range 1 ^ ^ ^ 10"^. Below, we show that a similar trend is 
found for the magnetic support. 

Figs. [9] and [To] show the dependence of the magnetic 
field amplification on the divergence and the shear, where 
velocity derivatives are multiplied by the grid scale A to 
level out the different refinement levels. Let us first con- 
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Figure 8. Plots of the positive (solid lines) and negative (dashed lines) components of the shear-induced and compressive magnetic field 
amplification (see Eq. [sj as functions of overdensity at time t = 0.1% (light colour) and t = O.Stff (full color). 




Figure 9. Plots of the magnetic field amplification by contraction as functions of negative divergence at time t = 0.1% (light colour) 
and t = O.Stff (full color), as in Fig.js] 
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Figure 10. Plots of the positive (solid lines) and negative (dashed lines) components of dynamo term in Eq. ([s} as functions of trace- free 
rate-of-strain scalar at time t = O.ltff (light colour) and t = 0.5% (full color), as in Fig. [s] 
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Figure 11. Averages of A-^pAthermi as functions of the thermal 
pressure (top), A-^pAmagn ± as functions of the magnetic pressure 
(middle), and A-^pAturbi as functions of the enstrophy densities 
(bottom) for a weak magnetic field {(3q = 20) at t = O.Stff. As 
in Fig. |2] positive and negative components are shown as solid 
and dashed lines, respectively, the light colored lines are obtained 
form t = 0.1%, and identity functions are indicated by dot-dashed 
lines. 



sider the weak-field case (/3o = 20). In the case of compres- 
sive amplification, only the contribution from negative di- 
vergence is plotted. Except for strongly negative values of d 
at time t = 0.5tff, where the effect of gravitational compres- 
sion becomes manifest, {—B'^d)d is nearly proportional to d. 
This implies that regions of higher divergence are not associ- 
ated with systematically higher magnetic fields. In a similar 
way, {—BiBjSij^)\s*\ varies roughly linearly with the rate- 



0.001 0.1 10 1000 10^ 10^ 

2 

Figure 12. Same functions as in Fig. |11| for a strong magnetic 
field (/3o = 0.2). 



of-strain scalar \S*\. Thus, the magnetic field cannot remain 
attached to turbulent structures that efficiently amplify the 
field. In other words, there is no strong backreaction of the 
field on these structures. This is not a odds to the frozen-in- 
motion in ideal MHD because neither shear nor divergence 
are properties carried by fluid elements. On the other hand, 
the field tends to be captured in the collapsing gas and this 
is why we see the super- linear increase for —dA > 1 af- 
ter 0.5 free-fall time scales. For the strong field (/3o = 0.2), 
the behaviour of the compressive amplification is qualita- 
tively similar, but the dynamo action is markedly different. 
The graph of {—BiBjS*]^)\s*\ is significantly flatter than for 
/3o = 20, but the amplification rises steeply for strong shear. 
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In contrast to fio = 0.2, this indicates that back-reaction 
effects cause a strongly non-Unear interaction between the 
turbulent gas flow and the magnetic field. There is only lit- 
tle net amplification though, as the positive and negative 
components nearly balance each other. 

As motivated in Sect. [4] the mean thermal, magnetic, 
and turbulent support functions are plotted as functions of 
P, Pm and lA^pcj^, respectively. The results for /3o = 20 
are shown in Fig. [TT] At t = 0.5tff, the thermal support 
(top panel) is qualitatively similar to the case of hydrody- 
namic turbulence fsee Fig.|2|, but the range of pressures, for 
which the mean A^pAtherm+ approximately matches P, is 
very narrow. This is probably due to the significantly lower 
resolution limit. The magnetic support function (middle 
panel) shows interesting properties. The overall behaviour 
is quite similar to the thermal support, which suggests that 
the magnetic support is basically pressure- like, althogh the 
non-diagonal Maxwell stresses are included in the support 



function defined by Eq. (16). This becomes also apparent in 



the phase plots of the thermal and magnetic support func- 



tions (Figs. A5 and A6). Only the negative magnetic sup- 
port (compression due Maxwell stresses) is relatively large 
at lower magnetic pressure. For high magnetic pressures, 
Pm ^ 1000, the magnetic field produces strong support, 
which becomes comparable to the thermal support. Thus, 
magnetic fields have a significant influence on the stabi- 
lization of the gas against gravity even if the field is rel- 
atively weak. Negative turbulent support (bottom panel) 
again dominates, although the linear relation found for hy- 
drodynamic turbulence is not as closely matched in the mag- 
netic case. For the strong field (/3o 



12 



shows 



0.2), Fig. 

qualitatively similar results as for the weak field, except for 
significantly less magnetic support. This result seems at first 
counter-intuitive, but it can be understood by considering 
the probability density functions of the magnetic field (see 
Fig. 13 in Collins et al. 2012j). Although the mean field is 
stronger for /So — 0.2, it turns out that the magnetic field 
fluctuations have a much wider tail toward high field inten- 
sities for /3o = 20. This is why the local support function, 
which is mainly caused by strong fluctuations, tends to be 
stronger in the latter case. 

Profiles of A/47vGpo6 against the overdensity 1 + ^ are 
plotted in Fig.|13| At the early instant t = O.ltff , the support 
gradually decreases with the density and approaches unity 
for S above 200. Clearly, supersonic compression is stronger 
than the magnetic support of the gas, although the mag- 



netic compression rate An 



exceeds At, 



at moderate 



overdensities in the case /3o = 0.2. In the advanced stage of 
gravitational collapse (t = 0.5tff), there is large positive sup- 
port for overdensities below 10^. K/AtiGpq5 becomes smaller 
than unity only for ^ - 10^. The main support is thermal, 
particularly in the strong-field case (/3o = 0.2). This corre- 
sponds well to the trends for the pressure-like quantities in 
Figs. [TT] and [12] For /3o = 20, however, the magnetic support 
is almost as strong as the support by thermal pressure In 
both cases, the positive component of the turbulent support 
is much smaller. The negative turbulent component tends 
to be comparable to the negative magnetic component, al- 
though the former dominates for high densities {8 > 10^) 
at time t = 0.5tff. The large impact of magnetic fields on 
the support of the gas for f3o — 20 is further demonstrated 
by the phase plots shown in Figs. |A7| to |A9| The trend of 



Atherm +/47rG'po^ to bccomc smaller than unity is particu- 
larly pronounced in this case. Apart from that, the distri- 
butions of the thermal and turbulent support functions vs. 
density and the turbulent support spectra (not included in 
this paper) are qualitatively similar to the hydrodynamical 
case. 

The profiles computed from the root-grid representation 
of the data are plotted in the bottom panels of Fig. |13| In 
this case, the support drops below unity for densities above 
the critical density, which is set by the resolution criterion 



based on the local Jeans length (see Collins et al. 2012). 
As discussed in Sect. [4] this indicates that the high-density 
gas is dominated by gravity when smoothed out over the 
length scale of the root grid. In contrast, the pressure of 
smaller structures at higher refinement levels can locally op- 
pose gravity up to very high densities, but these structures 
are not sufficiently resolved relative to the Jeans length in 
these MHD models. A further important result is that mag- 
netic support is relatively weak in the range of densities, 
where A/AnGpoS < 1. Consequently, the relatively strong 
magnetic support that can be seen for the full AMR data in 
the case fio = 20 stems from local field compression and fold- 
ing, which largely average out over the root grid cells. On 
the other hand, the positive contribution from turbulence 
becomes comparable to the thermal support at the highest 
densities and exceeds the support by magnetic fields. Thus, 
the statistics of Aturb/47rG'po plotted in Fig. 13 show that 



vortices on the root-grid can play an important role in pro- 
viding support against gravity in the most dense gas, while 
shock compression generally dominates on the smaller scales 
of the higher refinement levels. 

Since gravity breaks the scale invariance of turbulence, 
the properties of collapsing turbulent structures in a numer- 
ical simulation are inevitably resolution-dependent. If the 
grid scale decreases, self- gravitating overdense gas will col- 
lapse to higher densities, unless the physical length scale 
on which the collapse ends by whatever mechanism is re- 
solved. In the context of star formation in the interstellar 
medium, the end of the collapse would be reached when the 
gas becomes optically thick. Since this is beyond the achiev- 
able maximum resolution for a box size of the order of par- 
sees, the collapse could be stopped numerically by dump- 
ing mass into sink particles beyond a critical density that 
is set by the resolution limit. For the reasons outlined in 
Sect. [T] however, no sink particles were applied in our sim- 
ulations. Having said that, gradients of the pressure, the 
velocity, and the magnetic field become steeper as smaller 
length scales are resolved, which give rise to stronger fluc- 
tuations and thus enhance the magnitude of the local sup- 
port. It is difficult to see a priori how this affects the ratio 
A/47vGpoS. So the question arises how robust the trends we 
observe are if the resolution is varied. We investigate this 
question by changing the resolution of the root grid, while 
keeping the number of refinement levels constant. Then the 
size of the densest structures that can be resolved at the 
highest refinement level decrease proportional to the root- 
grid resolution. Apart from that, the overall turbulent ve- 
locity field becomes coarser, corresponding to a higher nu- 
merical viscosity. We plot the results for A±/47vGpoS for the 
three different root-grid resolutions 128^, 256^, and 512^ in 
Figs. ^] and |15| The Truelove number is 16 for all simula- 
tions. Clearly, there are systematic differences in the positive 
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Figure 13. Plots of the positive (thick solid) and negative (thick dashed) components of the total support relative to the gravitational 
compression rate in Eq. \12\ as functions of overdensity. The left and right column of plots show the cases of weak and strong magnetic 
fields, respectively. The plots on the top are obtained for t = 0.1%. For t = 0.5tff, the results from both the full AMR data (middle) 
and the root-grid data (bottom) are plotted. In each plot, the dominant components of the thermal and turbulent support functions and 
both components of the magnetic support are indicated by thin lines. 



support at high densities. Structures in refined regions re- 
sist gravity up to higher densities for increasing resolution. 
This is expected, according to the discussion above. Owing 
to the randomization and the strong intermittency of lo- 
cal turbulent structures, the spikes in the support functions 
differ substantially. Nevertheless, the trend with density is 



qualitatively similar for the three simulations. In this regard, 
time-averaging would be advantageous, but this is not fea- 
sible because self-gravitating turbulence is not statistically 
stationary. If we consider the support computed from the 
smoothed fields on the root grid, we also see a systematic 
drift of the density range with A/AnGpod < 1. This implies 
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Figure 14. Plots of the total support as in Fig. [13] for different 
resolutions of the root-grid in the case of weak field (/5o = 20) and 
t = 0.5%. The statistics is plotted both for the full AMR data 
(top) and the root-grid representation of the data (bottom). 



that the structures captured by the refined grids at higher 
resolution indeed have a noticeable influence even when they 
are smoothed out over the root-grid scale. The increase of 
the support with resolution at a given overdensity in the 
range 10 < 5 < 10^ can be attributed to the steepening of 
the gradients in Eq. (13), as the effective numerical viscos- 



ity becomes less. Naturally, for the reasons outlined above, 
convergence cannot be achieved as long as the gas is not 
prevented from collapse to ever higher densities. However, 
by conservative down-sampling of the supp ort functions for 
the higher resolutions to 128^ (see Fig.nGj), we also see that 
the main difference can be attributed to the averaging from 
larger to smaller scales. Thus, the general trends of the sup- 
port we outlined above are not substantially altered by the 
numerical resolution. 
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Figure 15. The same plots as in Fig.fMlfor /3o = 0.2. 
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Figure 16. Plots of the root-grid representation of the total sup- 
port down-sampled to 128^ for /3o = 0.2 and t = O.Stff . 



6 CONCLUSIONS 

In the highly non-linear, turbulent regime, we propose to in- 
fer the local contributions of turbulence, thermal pressure, 
and magnetic fields to the support of the gas against grav- 
ity from the dynamical equation for the rate of compres- 
sion, —Dd/Dt (Schmidt 2009| . We performed a statistical 
analysis of the source terms, i. e., the local support A and 
the gravitational compression rate 4:7TGpoS, where po is the 
mean density and S the relative density fluctuation, for three 



AMR simulations of self- gravitating turbulence. One simu- 
lation is purely hydrodynamical ( Kritsuk et al.pOlla ), the 
other two are MHD simulations with different mean field 
strengths ( jCollins et al. 12012 ^. In the following, we summa- 
rize the main results of our study: 



(i) The statistics of the local support function (14) does 
not provide evidence of positive turbulent support. On the 
contrary, the net turbulent support is negative at all den- 
sities and for arbitrarily high vorticity. This result is ob- 
tained both for hydrodynamic turbulence with mixed fore- 
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ing and MHD turbulence with purely solenoidal forcing at 
high Mach numbers. The spectra of this function also show 
a dominant negative component for all wave numbers. The 
negative component can be interpreted as negative pressure 
produced by the strain of the turbulent velocity field. For 
hydrodynamical turbulence, we find a simple linear relation 



significantly modulate the volume integrals of the thermal, 
turbulent, and magnetic energy densities (see, for instance, 
Lequeux]|2005 ), depending on conditions that vary in space 



(Eq. 20) between the vorticity and the negative pressure. 



Magnetic fields cause deviations from this linear relation. 

(ii) The largest contribution to the support of the gas 
comes from fluctuations of the thermal pressure, also in the 
case of MHD. For isothermal gas, this means that compres- 
sion rate is mainly given by the density structure of the 
gas because both the thermal support function ( 15 ) and the 
gravitational compression rate 47vG{p — po) are solely deter- 
mined by p. 

(iii) Magnetic fields contribute significantly to the sup- 
port. For an initial field with /So = 20 (corresponding to 
P ^ 0.2 for statistically stationary turbulence), the positive 
component of the local magnetic support function (16) is 



comparable to the thermal support. However, the support 
does not further increase for a stronger field with /3o = 0.2 
~ 0.03). The reason for this is the wider tail of the PDF of 
B for larger fio . While turbulence causes the initial compres- 
sion of the gas, the turbulent support function in overdense 
gas is small relative to the magnetic support function after 
0.5 free-fall time scales (except for very high, underresolved 
densities) . 

(iv) Gravitational collapse locally enhances the magnetic 
field strength. In addition to the compression by gravity 
and shocks, a roughly equal amplification is caused by shear 
(small-scale turbulent dynamo). Consistent with our find- 
ings for the magnetic support, the amplification in high- 
density gas tends to be stronger for a relatively weak mean 
field (/3o = 20). 



(v) The local support function, as defined by Eq. (13), 
varies between values greater or smaller than the gravita- 
tional compression rate. For moderate overdensities, we find 
a negative net support, i. e., gravity and supersonic turbu- 
lence cause the gas to contract. At higher densities, however, 
positive pressure support becomes strong and the average 
ratio of the support function to the gravitational compres- 
sion rate is much greater than unity, with strong fluctuations 
due to intermittent structures. Gravitational compression is 
dominant only for very high densities at the resolution limit. 

(vi) We compared our results to smoothed data, which 
can be obtained from the root-grid representation. The main 
effects of the smoothing are that the total support becomes 
small compared to gravitational compression for much lower 
densities, marginal turbulent support is found at the high- 
est densities, and the support by magnetic fields is reduced. 
Thus, it appears that the gas is self-gravitating when av- 
eraged over large volume elements, while smaller structures 
that are only resolved at higher refinement levels are sup- 
ported by strong fluctuations of the pressure and the mag- 
netic field. 



The behaviour that emerges from our statistical anal- 
ysis is unexpected in several aspects, most noticeably, the 
lack of a positive turbulent support. Global support by tur- 
bulent pressure is clearly an implication of the generalized 
virial theorem for an isolated object. For a smaller subregion 
inside a turbulent cloud, however, the boundary terms will 



and time. The local support functions we use in this arti- 
cle are source terms of the partial differential equation ([5| 
for the divergence of the velocity, which can be directly re- 
lated to the mass density source of the gravitational po- 
tential in the Poisson equation. The positive and negative 
variations of the local support function imply that a partic- 
ular fluid element frequently experiences contraction phases 
with A/47vGpoS < 1 and then crosses regions with strong 
positive support, where A/47vGpo6 > 1. Only if it accumu- 
lates enough density so that gravitational compression is 
dominant over a sufficiently long period of time, it will be 
torn into irreversible collapse. To improve our understand- 
ing of the mechanism of core formation, Lagrangian statis- 
tics might turn out to be valuable. The role of supersonic 
turbulence is ambiguous in this regard. Shocks can compress 
the gas to sufficiently high density, but also disrupt dense 
cores. On the average, the compression effect is dominant. 
The analysis of simulations of self-gravitating turbulence by 



Federrath & Klessen (2012) also suggest that the main role 



of turbulence is to enhance gravitational collapse. It is pos- 
sible though that non-turbulent, rotation-free flow produced 
by the pull of gravity biases the turbulent support toward 
compression because the negative divergence of such flows 
contributes to the rate-of-strain scalar. Currently, it is un- 
clear how to separate this from the genuinely turbulent ve- 
locity fluctuations. But the contribution of gravity to the 
turbulent support function appears to be small compared to 
the impact of shocks for several reasons. Firstly, we see the 
dominance of negative turbulent support already at early 
stages, where the gas only begins to contract under its self- 
gravity. Secondly, if the negative component of the turbu- 
lent support were strongly influenced or even dominated by 
gravity, it should decrease relative to the positive component 
towards low densities, as gravity becomes increasingly weak. 
This is also not the case. Moreover, one can clearly see in 
Fig.^that strongly negative turbulent support is associated 
with shocks. Thirdly, the velocity power spectra are not sig- 
nificantly affected by gravity ( Collins et al.|2012| Federrath 
k Klessen|[20T3t . 



In the case of hydrodynamical self-gravitating turbu- 
lence, the extremely high dynamical range of the simulation 
allowed us to identify a nearly self-similar regime with a neg- 
ative net support of magnitude |A| ~ AnGpod. This regime 
could be roughly described by the free-fall collapse solution 
of|Penstoii|( |1969| . However, we find a subdominant, yet non- 
negligible positive component of the thermal support. The 
idea of a hierarchy of quasi-virialized structures proposed 
by Biglari & Diamond (1988), on the other hand, does not 



appear to be consistent with our simulations, because this 
would imply a positive net support of the order of the grav- 
itational compression rate. 

We attempted to bridge the gap between local and 
large-scale support by computing the support functions from 
the smoothed data on the root grid of our AMR simulations. 
Since the refinement follows the gas density, this is roughly 
equivalent to integrations over larger volumes of overdense 
gas and, indeed, the positive component of the turbulent 
support tends to become larger. On the other hand, the 
statistics computed from the fine grid data reflects the dy- 
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namics of the collapsing cores, at least up to densities, for 
which the physics of self-gravitating turbulence is sufficiently 
resolved. It is an interesting question whether the different 
behavior of the support we see on the coarse root grid and 
on the refined grids could be carried over to observations 
of molecular clouds. Indeed, the root-grid resolution, which 
is about 0.01 pc roughly corresponds the pixel size of cur- 
rent CO maps (e. g.. Goldsmith et al. 2008). Some observed 
prestellar molecular cores have Bonnor-Ebert-like profiles 
with central densities up to several 10^ cm~^ (for a recent 



review, see Bergin & Tafalla 2007). For the corresponding 
overdensities ^ ~ 10 at the root grid of our hydrodynamical 
simulation, we find a positive net support A that is compara- 
ble to the gravitational compression rate 47vGpo6. This could 
mimic Bonnor-Ebert-like cores in quasi-hydrostatic equilib- 
rium. In stark contrast, however, the small-scale structure 
that is resolved by AMR reveals collapsing cores resulting 
from supersonic compression, while the thermal pressure 
support is subdominant. These considerations are supported 
by calculations of the instantaneous virial balance of clumps 
and cores, including the surface terms, by Dib et al. ( 2007| ). 
They show that clumps and cores are non-equilibrium struc- 
tures because surface terms are generally of the same or- 
der as the volume terms. In particular, they find structures 
that are undergoing supersonic compression, but not all of 
them are gravitationally bound. As a future extension of 
our study, clump finders or, alternatively, the dendrogram 



analysis technique developed by Rosolowsky et al. ( 2008| 
could be applied to constrain the statistics of the support 
functions to clump-like objects on different length scales. 
This might also help to interpret observational data in the 
light of high-resolution numerical data and to clarify which 
biases are introduced by relatively coarse observational reso- 
lutions. The dendrogram analysis of numerical and observa- 
tional data presented in [Rosolowsky et aT] ([2008 ) , for exam- 
ple, reveals scale-dependent discrepancies in the distribution 
of self-gravitating objects that are characterized by the virial 
parameter. 

Magnetic fields locally have a significant stabilizing ef- 
fect fsee also [Federrath Kless'eii||2012t , but we find that 
magnetic support does not dominate for large mean mag- 
netic pressure (/3 ^ 1). The stronger support in the case of 
a weak field (/3o = 20) can be explained by the turbulent am- 
plification in the initial turbulence simulation and the result- 
ing wider tails of the magnetic field strength in comparison 
to Po — 0.2 (Collins et al. 2012). In conclusion, the mean 



magnetic pressure is not a suitable quantity to infer the 
importance of magnetic fields in stabilizing self- gravitating 
gas. By computing the magnetic support function for the 
smoothed root-grid data, we find a weaker effect in com- 
parison to both the thermal and turbulent support func- 
tions. This suggests that the magnetic field is significantly 
amplified in local, collapsing structures. Indeed, we find a 
large amplification effect by gravitational and supersonic 
compression at high densities, which is roughly balanced 
by shear-induced amplification. This is different from the 
results of |Sur et al]|2010 ) and [Turk et~al] ( |2012 ) , who run 
simulations of collapsing halos, where subsonic turbulence is 
solely produced by the gravitational collapse. Since turbu- 
lence is vortex- dominated on length scales below the Jeans 
length, the amplification by shear becomes strong compared 



to compressive amplification in this case if the numerical 
resolution is high enough. 

The importance of magnetic fields in molecular clouds 
has been notoriously difficult to pin down (see Crutcher] 
2012). There are observations that show fields are possi- 



bly weak, such as the column density power spectra study 
of [Padoa n et al.| ( |2004 | , as well as observations that they 



are dynamically important ( jLi et al.||2009| . Kritsuk et al 



( ^2011bJ argue that the super- Alfvenic nature of turbulence 
in molecular clouds is a natural outcome of the cloud forma- 
tion process. Our simulations show that the relative impor- 
tance of magnetic fields is dependent on density scale as well 
as averaging scale. Thus it is important to carefully consider 
both effects when interpreting observations using numerical 
data. 

In our resolution study, we demonstrate that the effect 
of numerical viscosity has a potential impact on the statis- 
tics of the local support function. However, a direct estimate 
of the effect of numerical viscosity, for example, as proposed 
by Pan et aT] (f2009 ) , is infeasible because the mean numeri- 
cal viscosity can only be estimated from ensemble averages. 
These can be approximated by global averages for uniform- 
grid simulations, but not for the different refinement levels 
in AMR simulations. An improvement could be made, how- 
ever, by computing the effect of unresolved turbulence with 
a subgrid scale model. Unfortunately, a subgrid scale model 
for highly compressible MHD turbulence is not available yet. 
For hydrodynamical turbulence, on the other hand, the sub- 



grid scale model by Schmidt & Federrath ( 201l| can be ap- 
plied. By adding the subgrid scale stresses to the equation 
for the compression rate, not only the effect of turbulent 
viscosity, but also the contribution form the turbulent pres- 
sure due to velocity fluctuations below the grid scale can be 
calculated. lapichino et al. (12011") show that this contribu- 
tion is small, yet non-negligible in numerical simulations of 
turbulence in the intergalactic medium. For supersonic tur- 
bulence in the interstellar medium, however, the unresolved 
fraction of the turbulent pressure can become comparable 
to the thermal pressure. So far, it is not clear whether this 
effect would enhance or further decrease the support of the 
gas. A further implication of our study is that the analy- 



sis of turbulent pressure support by jZhu et al. (2010) and 
[lapichino et al.| ( [2011[ ) falls short of the possibly dominating 
negative component, at least for the highly compressible tur- 
bulence in the intergalactic medium. It is clearly necessary 
to consider both the positive and the negative contributions 
to the local support. 

Our findings also bear consequences on theoretical and 
numerical approaches to the clump or core mass function 
in star-forming clouds. For example, the analytic theory 



by Hennebelle & Chabrier (2008) assumes that turbulence 
has a scale-dependent stabilizing effect, which follows from 
the incorporation of a power law for the turbulent pres- 
sure into the Jeans criterion. As shown by [Schmidt et al.[ 
( 2010| , this has a significant impact on the mass spec- 
trum of gravitationally unstable clumps in hydrodynamical 
turbulence simulations (without explicit treatment of self- 
gravity). In a similar way, an effective equation of state is 
used in the statistical excursion-set model for fragmentation 
in self- gravitating turbulent media, which was recently put 



forward by Hopkins (2012). From our discussion, however. 



it follows that an effective local Jeans mass for the clumpy 
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substructure of supersonic turbulence is difficult to consoli- 
date with the properties of the turbulent support function. 
This does not mean that the Chandrasekhar |1951 ) rela- 



tion for the effective speed of sound is invalid per se. Such a 
relation can be rigorously derived for hydrodynamical tur- 
bulence by means of scale decomposition (see Schmidt 
Federrath|2011 ) and, on length scales larger than the forcing 
scale, also for self-gravitating turbulence (Bonazzola et al. 



1992). But the fully non-linear dynamics of gas compres- 
sion, which is determined by the local balance between ther- 
mal, turbulent, and magnetic support versus gravity, de- 
pends on the Laplacians and gradients of the pressure com- 
ponents. As opposed to the thermal and magnetic pressures, 
our statistical analysis implies that large turbulent pressure 
does not give rise to a strong mean support against gravity 
for highly resolved self- gravitating turbulence. As a conse- 
quence, the Lagrange identity form of the virial theorem. 



2(^kin + ^int) + 

substructure in molecular clouds. 
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the mean value of Y for X = Xi is given by 

{Y)i = ^ n dv{Xi, n) (j2 '^(^^' '^'^)) = 

k \ k / 

= ^(n+-n_)di>(Xi,n) /'^di>(Xi,n)'j = 

k \ k J 

= - . 

Note that {Y±)i / Y±- because the sums in the de- 
nominators are generally different. Analogous expressions 
hold for mass- weighted averages with dm{X,Y) in place of 
dV{X,Y). 



APPENDIX A: PHASE PLOTS 

In the following, two-dimensional distributions are shown, 
where the volume or mass fractions occupied by each bin 
are indicated by colour scales. In each case, the positive and 
negative contributions are plotted separately so that we can 
use logarithmic scaling in the plots (we do not show, how- 
ever, distributions of logarithmic variables). The distribu- 
tions are computed from the full AMR data sets at time 
t = 0.43tff for the hydrodynamical simulation and t = 0.5tff 
for the MHD simulations. The curves show averages of the 
quantity on the ^/-axis vs. the quantity on the x-axis normal- 
ized by the integrated volume or mass of those cells in which 
the quantity on the ^/-axis is non-zero. This normalization 
avoids a bias of the averages relative to the two-dimensional 
distributions. Even so, cuts through the distributions in y- 
direction can have very wide, asymmetric and irregular tails, 
which lead to pronounced fluctuations in the mean values. 
In mathematical terms, the distribution dV{X, Y) is a mea- 
sure of the differential volume occupied by particular values 
of the random variables X and Y. In the phase plots, we 
show the piecewise constant functions dV{Xi,Yk ±), which 
approximate dV{X, Y±). The mean value of Y±, given a par- 
ticular value X = Xi, is 

k \ k 

where it is understood that dV{Xi, Yk±) > only if Yfc ± > 
because, for example, all bins for which Yk < are unoc- 
cupied by Yfe + . 

In Figs. [2[ ]16| on the other hand, we normalize by the 
unconstrained volume or mass of each bin so that (A) = 
(A+) — (A_) is satisfied (see Sect.[4| and the profiles of neg- 
ative and positive components can be directly compared (the 
bin size is also smaller by a factor of 2 in these plots and, 
for as smoother representation of the data, curves are com- 
puted with Mathematica's higher-order interpolation func- 
tion rather than the linear interpolation between data points 
that is used by yt). This means that the decomposition of 
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Figure Al. Phase plots of of the normahzed thermal support against the thermal pressure (no magnetic field). 




Figure A2. Phase plots of normalized turbulent support against the enstrophy density (no magnetic field). 
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Figure A3. Phase plots of the thermal support relative to the gravitational compression rate (no magnetic field). 



Local support against gravity in magneto -turbulent fluids 



10' 




l.OOe-03 

l.OOe-04 

l.OOe-05 

l.OOe-06 

l.OOe-07 

l.OOe-08 - 

l.OOe-09 

Il.OOe-10 
l.OOe-11 
1.006-12 
l.OOe-13 



Figure A4. Phase plots of the turbulent support relative to the gravitational compression rate (no magnetic field). 





Figure A5. Phase plots of normalized thermal support against the thermal pressure (/3o = 20.0). 
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Figure A6. Phase plots of normalized magnetic support against the magnetic pressure (/3o = 20.0). 
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Figure A7. Phase plots of the thermal support relative to the gravitational compression rate (/3o = 20.0). 
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Figure A8. Phase plots of the turbulent support relative to the gravitational compression rate (/3o = 20.0). 
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Figure A9. Phase plots of the magnetic support relative to the gravitational compression rate (/3o = 20.0). 



